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Abstract 

The system of 4 differential equations in the external invariant satisfied by the 
4 master integrals of the general massive 2-loop sunrise self-mass diagram is solved 
by the Runge-Kutta method in the complex plane. The method, whose features are 
discussed in details, offers a reliable and robust approach to the direct and precise 
numerical evaluation of Feynman graph integrals. 
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1 Introduction. 



High precision measurements in high energy physics (or more in general in the determi- 
nation of particle properties) require more and more precise calculations of multi-loop 
Feynman diagrams to have sufficiently precise theoretical predictions to compare with. 

The nowadays widely accepted procedure of expressing radiative correction calcula- 
tions in terms of a limited number of master integrals (MI) [|IJ reduces the problem to 
the careful determination of these quantities. The method has also the advantage that, 
with a correct bookkeeping of the recurrence relations arising from integration by parts 
identities, the MI of a given problem can be reused in more complicated calculations. 

The analytical calculation of MI, in terms of the usual polylogarithms and their gener- 
alizations, is in general possible only when the number of different scales (internal masses 
and external momenta or Mandelstam variables) is small, like in QCD calculations, where 
all masses are set to zero or in the QED cases, where only the electron mass is different 
from zero, or when the external variables are fixed to particular values (zero or mass shell 
condition). Another possibility of big help in analytic calculations is sometimes offered 
by the exploitation of particular simplifying conditions, like the smallness of some ratios 
of the parameters allowing the corresponding expansion. 

In the general massive case, relevant in the electroweak theory, the number of param- 
eters prevents from obtaining results in the usual analytic form already in the case of the 
2-loop sunrise self-mass diagram shown in Fig.l. 




Figure 1: The general massive 2-loop sunrise self- mass diagram. 

This diagram has indeed a long history of investigation and its MI were even rec- 
ognized to be expressible in closed form as a combination of four Lauricella functions, 
a special class of generalized hypergeometric series (and earlier references therein). 
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The method provides efficient multiple series expansions for the regions of small |p 2 |, i.e. 
\p 2 \ < max(m 2 ), and of large \p 2 \, i.e. \p 2 \ > (mi + 1712 + m^) 2 , but some problems arise 
in the intermediate region. 

Great efforts were therefore devoted to investigate the properties in the special points 
(i.e. p 2 = 0,oo, pseudothresholds and threshold). The analytical expansions of the MI 
at and 00 are given in and ||; the values at pseudothresholds and threshold in [f|; 
the analytical expansions at pseudothresholds are in ||; a semi-analytical expansion at 
threshold is in || and also in configuration space technique in while the complete 
analytical expansions at threshold are presented in ||. 

For numerical evaluation purposes, it is possible to cast the general massive self-mass 
diagram as a double integral representation and in the particular case of the sunrise 
diagram in a single integral representation ||, || (and earlier references therein). The 
configuration space technique is also exploited in the numerical approach , pll . In a re- 
cent approach rearrangements of the integrand, driven by the Bernstein-Tkachov theorem, 
are introduced to improve numerical convergence [ IT| . A different and interesting method 



is the use of the recurrence relations as difference equations to numerically evaluate the 



MI pj. 



In the present paper we exploit the numerical evaluation of the four MI related to the 



general massive 2- loop sunrise self- mass diagram [13], using the differential equations in 
p 2 , obtained in J| and the Runge-Kutta method [11J] to solve them for complex values 
of p 2 . The interest is not limited to provide with a fast routine precise numerical values 
for all the four MI in the general massive case and for all the values of p 2 , but extends 
to the investigation of the reliability of the method, as it can be easily extended to the 
numerical evaluation of other less studied diagrams. 

In Section 2 the master differential equations are recalled and the analytic properties 
of the MI are reviewed. Section 3 contains a description of the method used to solve the 
system of differential equations and to determine the accuracy. In Section 4 the control 
tests and the comparisons with other values reported in the literature are presented. 
Finally in Section 5 our conclusions on the application of the method to present and 
further work are presented. 
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2 Analytical properties and behaviours of the MI. 



We use here the following definition of the four MI related to the general massive 2-loop 
sunrise self-mass diagram in n continuous dimensions and with fully Euclidean variables 



m i l 2 2 2\ 

F j (n,m 1 ,m 2 ,m 3 ,p ) 



8-2n 



((2tt)"- 2 ) 2 

d n h f d n k 2 — 2 

J (k\ + m\ ) Q1 



0')(jfe| + m 2 2 ) a ^\{p -h- k 2 ) 2 + m 2 3 ) a ^ 



,j = 0,1,2^1) 



and z = 1, 2, 3; a;(0) = 1, for j = 0; a^j) = 1, for j ^ i; a^j) = 2, for j = i. 

Wherever necessary to avoid ambiguities, the usual imaginary displacements m 2 — > 
m 2 — ie, where e is an infinitesimal positive number, are understood. 

At variance from ||, were the mass scale was given the value \i — 1, here we choose 

/i = m 1 + m 2 + m 3 , (2) 

which comes out to be the appropriate mass scale parameter for the numerical discussion. 
The expansion of the MI around n = 4 has the form || 



Fj(n, wii, m 2 ,m 3 ,p ) = C (n) 



i(-2) / 2 2 2 2\ 

j K,m 2 ,m 3 ,p) 



+ (^4) ^'^ ( m ?' ™ 2 ' m ^ P 2 ) + ^ (0) (ml ml mj, p 2 ) + 0(n - 4) | . (3) 



where the coefficient C(n) 



c( n ) = (2^) M r(3 



(4) 



not expanded, can be replaced by its value C(4) = 1, at n = 4, when multiplying a 
function regular in (n — 4). The coefficients of the poles in (n — 4) of F (n, ml m 2 , m 2 , p 2 ) 
are known to be HI 



tti(-2)/ 2 2 2 2\ 

F V '{m 1 ,m 2l m 3 ,p ) 



rp(— 1)/ 2 2 2 2\ 

F V '{m 1: m 2 ,m 3 ,p ) 



lfp ^ / 2 2 2\ 

8 |T + 2^ 1 + m2 + m ^ 

2 2 2 

m 2 log(^) + m 2 . log(^) + m 2 logfe 



(5) 



while those for Fj(n, m\, m\, m|,p 2 ), % = 1,2,3 are 



t-,(-2)/ 2 2 2 2\ 
,(-l)/^2 ^,2 ^,2 „,2\ ! 1- - m 



^K.m^mt^) = -- + -l g(-±). (6) 

From now on we deal only with the finite parts of the MI (i.e. we subtract the 
n — 4 poles) and we do not write anymore, for short, the arguments of the functions 
Fj°\m 2 , m\, m\,p 2 ) = Fj , unless we need to refer explicitly to them. 

The differential equations satisfied by the finite part of the MI expansions, given in 
[[J, can be written as 

2 9 p (0) _ p (0) 2 p (0) 2 p (o) 2 p (o) T 
p 2 D(m 2 ,m 2 ,m 2 ,p 2 ) Aif> = £ M^-if > + T, , i = 1,2,3 (7) 



j=0 



where the explicit form of the functions T ,Tj (polynomials of p 2 and mf, and loga- 
rithms of m 2 / fx 2 ) and Mjj (polynomials of p 2 and m 2 ) can be found in ||. The function 
D{m\,m 2 ,m\,p 2 ) is defined by 

D{m\,m\,m\,p 2 ) = (p 2 - p 2 fo ) (p 2 - p 2 psl ) (p 2 - p 2 ps2 ) (p 2 - p 2 ps3 ) , (8) 

Pth = ~( m i +m 2 + m 3 ) 2 , 

P 2 p si = -(mi+m 2 -m 3 ) 2 , 

P 2 P s2 = ~( m i -m 2 + m 3 ) 2 , 

P 2 P s3 = -( m i - m 2 - m 3 ) 2 , 

and vanishes at the threshold p 2 th and at the three pseudothresholds p 2 sl , p 2 2 , p 2 s3 - Indeed 
the values of p 2 at which the coefficients of the derivatives in Eq. (0) vanish, together with 
p 2 = oo, are the singular points of the differential equations, which we will call special 
points because special care is required in the numerical computation of the MI. 

The differential equations Eq.(^) allow a class of solutions wider than just the functions 
Eq. ([[]), but the initial conditions at p 2 = imposed by Eq.fll]) identify uniquely the 
solutions (actually the regularity at p 2 = 0, even without the explicit knowledge of the 
functions in that point, is enough to fix the solutions). Once the initial condition at p 2 = 
is fixed, one finds that the pseudo-threshold values p 2 = p psl , p 2 = p 2 s2 , p 2 = p ps3 are also 
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regular points [||, while the threshold value p 2 = p 2 h is a branch point || (in agreement, 
of course, with standard textbook results Ul5| ). 

It is convenient to use reduced masses and reduced external invariant 

= m i 2 = P 2 /q\ 

m i,r — . . j Pr — / i i \9 i V*) 

mi + m 2 + r«3 (mi + m,2 + va^y 

together with a dimensionless version of Fq°\ defined by 

Fo,r = 7 : : \i i (10) 

{m 1 +m 2 + m 3 ) z 

as the other master integrals are already dimensionless, the values of all the functions 
are now pure numbers. In terms of the new variables p 2 ,mi^ r the threshold is located 
a t Pth,r = Pth( m i,r) = — 1 and the pseudo-thresholds are in p 2 slr = P ps i( m i,r) , Pp S 2,r = 
Pps2( m i,r),pl S 3 t r = Pps3( m i,r)- We also recall here that according to the Euclidean definition 
our p 2 is positive for space-like p, and negative for time-like p. 

Typical plots of the real and imaginary parts of the MI for two sets of masses are 
shown in Figures and |5|. The plots are obtained by means of the FORTRAN program 
described in the next two sections and each consists of 6000 points per function calculated 
with the relative precision of 10~ 6 . 

The behavior of the functions depends of course on the values of the masses, however 
some of their properties are quite general, as it appears from the following discussion. The 
real part of F %. as seen in figures |2| and ^| has one local minimum and one local maximum 
for finite p 2 and goes to ±oo for p 2 —* ±oo. We do not have analytically the exact position 
of the local extrema, but as they appear to lie outside the region — 1 < p 2 < 1 they can be 
found approximately with the help of the asymptotic expansion (large p 2 ) for Fq°J , which 
is I 

< = | (l0g(p?) - f ) + ^ log 2 ^) E <r ~ ^ E < l0g«) + " " " • (11) 

The asymptotic behavior of Re F Q ° is obvious from Eq.flTTD, while for the positions 
of the maximum and minimum one finds approximately from the first term only, which 
is also the leading, p 2 max = —9.5 and p 2 min = 9.5, independent from the mass values. 
The approximate values of the function at those points, again from the first term only, 
are Re FqJ (p 2 )max ) = 0.3 and Re F^} {p 2 min ) = —0.3, also independent from the mass 
values. Taking into account the asymptotic behavior one expects at least three zeros of 



5 



o.s 



o.s 



O 4 



O 2 









1 










/ i 












































(3) 












(oy/ 


















1 







-30 -20 -10 o 10 20 30 

Figure 2: Plots of Re Fq} (labeled as (0)) and Re f/ 0) (labeled as (i)) function of 
p 2 r for mi = 2, m2 = 1, = 4 and /i = ni\ + m2 + m^. 



the function, provided that in non asymptotic region there are no additional extrema 
(which is actually the case). The second derivative of Re Fq° at threshold is infinite ||, 
but it does not change sign at that point, even if the position of the flex point of Re Fq^ 
is not far from the threshold. 

The other MI go one into the other by the exchange of the values of the related masses. 
From the expression of their expansion for large p 2 r 

if } = ~ logU 2 ) + ^ log( Pr 2 ) (log(mJ r ) + l) + • • • , (12) 

we see that their real parts all go to —00 in both asymptotic regions pi — > ±00, however the 
position of the maximum cannot be obtained just from the first terms of the asymptotic 
expansion, as it is positioned in the region of small p 2 r . The analytic expansions at 
threshold || show that the derivatives of the Re F^ are infinite, but they do not change 
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Figure 3: Plots of Re F r (labeled as (0)) and Re F> (labeled as (i)) as a function of p 2 . 
for mi = 1, rri2 = 9, m-^ = 200 and /i = mi + m2 + m^. 



sign exactly at that point. 

The imaginary parts of all the functions plotted in figures f| and |5| exhibit no com- 
plicated structure and their asymptotic behaviors can be simply deduced from Eq.(p]) 
and Eq. (|T2"D . Observe that, due to our definition of p 2 ., the proper analytic continuation 
for time-like p, hence negative pi = —\pl\, is obtained by giving a positive infinitesimal 
imaginary part to — p 2 r1 so that \og(p 2 r ) — * log(p^ — ie) = log \p 2 r \ — in (at variance with the 
default option of FORTRAN compilers). 



3 The numerical method. 

For the numerical solution of the system of differential equations we use fourth-order 
Runge-Kutta method [Q. The method starts from the known values of the solutions in 
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Figure 4: Plots of Im (labeled as (0)) and Im (labeled as (z)) as a function of p 2 
for mi = 2, m2 = 1, m3 = 4 and /i = mi + m2 + m^. 



an initial point, then it calculates the values of the solutions in a nearby point at distance 
A with an expansion in A based on the differential equations, omitting terms of order 
A 5 . Repeating the procedure along a path of length L in N steps, so that the step is of 
length A = L/N, a relative error of approximately iVA 5 = L 5 /N 4 is accumulated, and the 
requested accuracy is obtained by a suitable choice of L and N. This method is known for 
its robustness, and indeed it works quite well in our case allowing us to obtain a relative 
accuracy of 10~ 10 — 10~ 12 (the FORTRAN program is written in double precision) within 
reasonable CPU time (see discussion at end of this section). More sophisticated methods 
exist and could be implemented, but the simplicity of the used one has the advantage of 
a better control on the accuracy. 

To obtain the four MI related to the general massive 2-loop sunrise self-mass diagram 
we use the system of four linear differential equations given in Eq.fl7|). For the necessary 
initial conditions we use the values of the MI in the special points, where the differential 
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Figure 5: Plots of Im (labeled as (0)) and Im (labeled as (z)) as a function of pi 
for mi = 1, m>2 = 9, = 200 and fi = ni\ + m2 + 



equations simplify, allowing the analytic calculation of the MI; but as also the coefficients 
of the derivatives vanish there, to start the numerical evaluation from that points the 
values of the first derivatives must be provided as well. We use for that purpose the 
analytic values presented in || ^ [8|. Starting from pi = 0, numerical instabilities arise 
when approaching any of the pseudo-thresholds; therefore to obtain the value of the MI in 
the proximity of the threshold p\ h r = — 1 or of the pseudo-thresholds Pp Sljr , p\ a 2, T -> Pps3,n we 
take the threshold or pseudo-thresholds themselves as the starting points of the numerical 
evaluation. Further, for very large values of pf, it is more convenient to start from the 
value x = 4r = 0, so that the asymptotic expansion at large pi is also needed. 

Pr 

We use the Runge-Kutta method in the complex plain of pi; the initial condition 
guarantees as in the real case the uniqueness of the solution, provided that we do not 

-1 



cross the cut [16|, which extends in the present case from the threshold (p1 hr 



to 



-oo along the real axis. As already remarked, due to use of Euclidean variables in Eq. ([!]), 
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the proper sign of the imaginary part of the solutions is obtained with a path laying in 
the lower half complex plane of p 2 . Using complex variable has the advantage 

that the initial conditions at p 2 = can be used to obtain the solutions everywhere, with 
a path which does not approach too much pseudo-thresholds and threshold. This feature 
is relevant, as the method reaches much faster the required accuracy, when starting from 
p 2 = 0, rather than from anyone of the other special points - even if the result is of course 
independent of the chosen path. 

The program is organized as an independent subroutine, whose arguments are the 
input values of the masses rrii, of p 2 (which is real) and for the required accuracies, and the 
output values of the MI with their errors. Actually the input accuracies refer separately 
to the real part, the imaginary part and the absolute value of the functions. However the 
accuracy of the imaginary part is controlled by the program only when not vanishing, i.e. 
for p 2 < —1, but its value is in any case an indication of the precision of the result. The 
program starts usually from the initial conditions at p 2 = 0. However, if > 900, it 
starts from initial conditions at x = 1/p 2 . = 0, while if \p 2 — p 2 s \ < 0.001, where p 2 is the 
threshold or any of the pseudo-thresholds, it starts from the initial conditions at p 2 . 

When starting from p 2 = 0, the fourth-order Runge-Kutta method is applied to the 
system of equations Eq.(|7]) and the required final value of p 2 is reached following the 
rectangular path in the complex plane: (0,0), (0,-0. l),(^,-0.1),(p^,0). When pointing to 
time-like values p 2 < 0, this path has the merit of avoiding the numerically troublesome 
points of the threshold and pseudo-thresholds. The same path is also used for reaching 
space-like values p 2 > 0, although no special points occur along on the real axis, as it 
turns out that the requested precision is usually reached much faster along a complex 
path than along the real one. 

The fourth-order Runge-Kutta for the system of equations Eq.(|7]) is also used when 
the starting point of p 2 is one of the pseudo-thresholds Pp Sltr , Pp S 2, r i Pps3,n along a similar 
rectangular path: (p 2 ps , n 0), (p 2 ps ^ -0.1), (p 2 , -0.1), (p 2 r , 0). 

When starting from p 2 = — 1 (the threshold) the same system of equations cannot 
be used, as the first derivatives of the master integrals Fi, i — 1, 2, 3 are infinite at that 
point. Instead, we introduce new functions FQ h r „ F\ h through the definitions 

F 0°r = V /m l,'' m 2/ m 3; X th l °g( x th ) + F^r 

= ^ VgEgvgE Xth log(xa) + , , = i, 2 , 8 , (13) 
lo m i} r 

where x t h = P 2 - + (mi if + m 2jr + m 3>r ) 2 = p 2 + 1, we generate algebraically the system of 



10 



differential equations which they satisfy, and then we solve numerically that new system 
within the program. We do not report here the new system of equations as one can 
easily obtain it from the Eq.(^) using Eq. (|T3|) . The function does not have the same 
problems of the fI°\ i = 1,2,3, but the subtraction is performed anyway to simplify 
the equations. The required value of pf, is then reached along the triangular path in the 
complex plane (-1,0), ( (p 2 r - l)/2 ,-0.01), ( p 2 r ,0). 

In the asymptotic region > 900 we perform the change of variables pi — > x = I /pi, 
then we subtract from the MI the terms not vanishing at x = (the original MI are indeed 
divergent at x — > ) and finally we write a system of equations for the subtracted MI 
Fft and F t as , i = 1, 2, 3 defined as 

4? = | (logtf) - 4) + 4 E < - 4 log(p?) E < tog(m? ir ) 



32 V ov " r/ i / 32 " " ' ^ " l() 
1 

32 



1 ^mJ r (5-61og(mJ r )+log 2 K r )) + F as 



r 



i=l 



^ (0) = -4log 2 (p r 2 ) + ^log(p r 2 )(logK r ) + l 



32 ° v " r/ 16 

+ 1 (-1 - 4 log(m 2 r ) + log 2 (m 2 r )) + F? . (14) 

Again the new system of equations for Fq s t , F°- s can be obtained in a simple way 
substituting Eq.flI4|) into Eq.(^). The numerical solution is then obtained in the variable 
x along the complex triangular path (0,0), ( x/2 ,-0.01),(x,0). 

The errors assigned to the final results of the Runge-Kutta method are estimated by 
comparing them to the results obtained with a number of steps 10 times smaller then for 
the final results. The difference of the two results is taken as the estimate of the absolute 
error. To account for the cumulated rounding error, we estimate the relative error in a iV 
step calculation as y/N x 10 -15 , (as the program works in double precision), and then take 
the cumulated rounding error as the relative error times the value of the result. We finally 
take the sum of the absolute error and the cumulated rounding error as an indication of 
the error in the result. 

The initial number of steps iVj for \p%\ < 1 is taken to be Ni = 2/min(accuracies) 
(where min(accuracies) is the smallest of the accuracies required in calling the routine) 
and Ni = 2\pl\/min(accuracies) for \p^\ > 1, but is set to Ni = 20 if the number comes out 
smaller then 20. If the required precision is not reached the number of steps is increased 
by a factor 4, the system is solved once more and the procedure for estimating the error 
is repeated. For high required accuracy it might happen that the estimated error grows 
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yl Re jff Im Ff) Re pf Im Ff_ 

-1000. -113.349786296(3) 96.922241476(2) -0.4811355157(1) 3.255217395(1) 
-30. -0.230629539580(2) 2.31604333072(1) 2.86686026427(1) 2.46075767262(1) 
-15. 0.044413679180(8) 0.964565210432(6) 3.47247970401(2) 2.23096593908(2) 
-1.5 -0.2536003902785(5) 0.0032295545873(1) 5.03180993338(1) 0.564979990606(8) 



Table 1: The benchmark values of F^} and F^ for masses mi = 1, m 2 = 9, m 3 = 200 
and /i = mi + m 2 + m 3 . 

pi Re F 2 (0) Im F 2 (0) Re F 3 (0) Im F 3 (0) 

-1000. -1.4689531571(1) 2.393163661(1) -0.81238274708(3) 1.17931522036(6) 
-30. 0.880149099634(4) 1.624323639603(8) 0.1097269441157(5) 0.496397637240(4) 
-15. 1.26943734837(1) 1.42094925579(1) 0.195266696001(4) 0.366097572324(2) 
-1.5 2.066516140875(4) 0.239828506599(7) 0.2185459121886(5) 0.0195952353591(3) 



Table 2: The benchmark values of F^ and F^ for masses mi = 1, m 2 = 9, m 3 = 200 
and /i = mi + m 2 + m 3 . 



when the number of steps is increased (because of an accumulation of the rounding errors, 
etc.). In that case the program gives out the best result (i.e. the one with the smallest 
error). It may also happen, in the case the accuracy obtained in a given step is almost 
equal to the required one, that in the next step the accuracy obtained is much higher then 
the required one. 

Typical running times on PC with Intel Pentium III (1GHz) CPU are the following: 
for required accuracy of 10~ 7 a fraction of a second for \pl\ ~ 0.2, 2 seconds for \pf\ ~ 2 
and 8 seconds for \p%\ — 30; for required accuracy of 10~ n 20 seconds for \p%\ — 0.2, 3.5 
minutes for \p^\ ~ 2 and 59 minutes for \pf \ — 30. 

The program is available from authors upon request and we report in Tables [l],|2| and 
|3| a few results, which can serve as a benchmark. The reported results were all obtained 
asking the accuracies to be 10~ n . In Table ^| the values at pi = 0, at threshold {pi = — 1) 
and at three pseudo-thresholds (pi = p 2 ps2 ^ r -, V%2,,r) are calculated from the known 

analytical results || f|, ||, §] incorporated into the program, therefore no error is indicated 
(the imaginary parts vanish for those values of pi). 



12 



Pr fV £l £2 £3 



-1 


-0.279454902855371 


4.83093725524177 


1.89253423642110 


0. 


.185424043556224 


-0.99 


-0.2798928396415(5) 


4.80520291678(1) 


1.884686031133(5) 


0. 


.1846128537367(4) 




-0.280281633048667 


4.78848044412341 


1.87869903011051 


0. 


.183940826371101 


-0.9 


-0.2836780811878(5) 


4.68669935479(1) 


1.836380342548(8) 


0. 


.1786103121374(3) 




-0.286233415451605 


4.62888926299134 


1.80974413459626 


0. 


.174898795219002 


-0.825 


-0.2866587055221(5) 


4.620058473554(8) 


1.805561498336(6) 


0. 


.1742958493924(4) 


.Ppsl,r 


-0.286906928933491 


4.61498769104262 


1.80314761916276 


0. 


.173945546345814 


-0.8 


-0.2876221116285(5) 


4.60069487253(1) 


1.796297920597(6) 


0. 


.1729424293568(4) 


-0.1 


-0.3101507246241(3) 


4.261426874519(4) 


1.619837556072(1) 


0. 


.1432933391084(2) 





-0.312816604092084 


4.22788075922252 


1.60134292154365 


0. 


.139821925866842 


1.0 


-0.3340476235037(6) 


3.970691522343(6) 


1.455289136414(2) 


0. 


.1103513887593(2) 


30.0 


0.323213333716(3) 


2.33632892937(2) 


0.425519391740(4) 


-0. 


178170159306(2) 


1000. 115.539777092(3) 


-0.8050335305(2) 


-1.7888035846(2) 


-1. 


11980894482(3) 



Table 3: The benchmark values of FqJ, , F% and F% for masses mi = 1, m 2 = 9, 
m 3 = 200 and \i = m 1 + m 2 + m 3 . 



4 Tests and comparisons. 

Several checks were done in the past (|, || |J to verify that the analytical expansions of 
the master integrals in the special points, used within the numerical program, satisfy the 
differential equations and agree with the results existing in the literature. 

A remarkable feature of the extension of the RK-method to the complex plane is 
that it provides some natural self-consistency checks of the algorithm implementation. 
Starting from a special point and moving to a chosen value of p 2 with different paths 
in the p 2 -complex-plane, the values obtained for the master integrals should agree inside 
the errors of the method discussed previously. One has however to remember that paths 
chosen in the upper and lower half-plane, respect to the real axis, give opposite sign to the 
imaginary part of the master integrals for time-like values of the external invariant above 
threshold, p 2 . < —1. An even more complete test is to reach the same value of p 2 starting 
from different special points, hence following different paths, and compare the values of 
the master integrals at p 2 obtained along the various paths. If the values coincide inside 
the assigned errors the consistency between the differential equations, the expansions in 
the special points used as initial values, the implementation of the RK-method and the 
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algorithm for estimating the errors are cross-tested in a rather effective way. 

We have performed several of the mentioned checks in the different regions of p 2 
obtaining the requested agreement. 

The only published precise numerical results for the general massive case (all different 
non-zero mass values) are presented in B |TI|] , in the form of a combination of the general 
massive case with massless cases, to cancel the pole singularities in (n— 4). In our notation 
that combination is 

rri I 2 2 2 2\ Trl T7(°)f 2 2 2 2\ 17(0)/ 2 n 2 2\ 

T 12 3N{P ,m x ,m 2 ,m z ) = -16[ + F„ w (m^ m 2 , m 3 ,p ) - F^'(m v 0, m 3 ,p ) 

- F (0) (0,mim2,p 2 ) + F (0) (0,0,m2,p 2 )] , (15) 

which has also the property of being independent of n; the overall factor (—16) accounts 
for the different definition of the master integral in Eq.(|l]) and our p 2 corresponds to (— p 2 ) 
in |f2j and to s in [|TT|]. 

To obtain the values for Fq°\o, 0, m 2 ,p 2 ) we use the analytic formula presented in [|J, 
while for the values of F^\m\,m\,m\,p 2 ), F^\m\,Q,m\,p 2 ) and Fq°\o, m 2 ,, m 2 ,p 2 ) we 
use the present program. Although the value zero for the masses is not allowed, we have 
checked that the limit can be in practice reached numerically. Comparing the results 
obtained for the mass values from 10~ 6 to 10 -9 , we can estimate the error coming from 
having a mass not exactly zero, by taking the difference between the results obtained 
with the two smallest values used for the mass to be set to zero. As the error due to zero 
mass limit is sometimes comparable with the error due to the RK-method, we sum the two 
errors for each of the considered functions. The final error of T^n is the sum of the errors 
assigned by the algorithm to each of the contributing functions in Eq.flTsl). The larger 
errors (or less efficiency in calculations) come from the zero mass contributions, for which 



an approach based entirely on analytical expressions is in preparation |T7| . Furthermore 
the choice of equal values for two or even all the three masses, reduces the number of 
the independent equations in the system of differential equations, generating potential 
numerical problems, although less serious than those for the zero mass. 

In the values for Ti 2 3Ar(p 2 , m\, m 2 , tttJ) are presented, for different sets of the masses 
7711,1712,1713, and for the two regions of small \p 2 \ < (mi + m<i + m^) 2 and large \p 2 \ > 
(nil + m2 + ^3) 2 i n their Table 1 and 2 respectively. 

We repeat in Table § for the same values of the masses and p 2 the results of Table 1 of 
[§] for the multiple series (first entry), pushed to a large number of terms in some cases, 
and our results (second entry). The results are in excellent agreement. 

Also in Table 7 of |TT| the values for the same combination T 12 3at (there called S c ) are 
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mi 777,2 "73 


P 2 T123N 


P 2 T\23N 


3 3 10 


-9 -7.3129877443 
-7.3129877442(26) 


9 -6.93244055931 
-6.93244055924(50) 


2 3 10 


-20 -4.1493850173 

-4.1493850171(18) 


20 -3.63591843327 

-3.63591843320(95) 


2 2 10 


-25 -2.3353847298 

-2.3353847298(14) 


25 -1.9428452190 

-1.9428452191(10) 


1 2 10 


-30 -0.8117674738 

-0.8117674738(15) 


30 -0.6306847352 

-0.6306847353(12) 


1 1 10 


-49 -0.3167501084 

-0.3167501085(24) 


49 -0.1950338472 

-0.1950338472(21) 


3 4 15 


-50 -7.9471022759 

-7.9471022760(33) 


50 -6.8270303849 

-6.8270303852(19) 


3 4 20 


-100 -6.0171476156 

-6.0171476159(59) 


100 -4.9485063889 

-4.9485063897(42) 


3 4 20 


-150 -6.3903568683 

-6.3903568686(73) 


150 -4.7506023184 

-4.7506023184(66) 


5 5 25 


-150 -14.5339444977 

-14.5339444982(87) 


150 -1.21816923644 

-1.21816923648(68) 


5 5 25 


-200 -15.0523063012 

-15.0523063010(95) 


200 -11.8790613597 

-11.8790613597(83) 



Table 4: Comparison for small \p 2 \ < (mi + ni2 + m^) 2 . In each box the first entry is 
the value of the multiple series of Table 1 of , the second entry is our result (the error 
in the last digits is enclosed in parenthesis). Our value of p 2 corresponds to — p 2 in 0. 
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p 2 


T\23N 


P 2 


T\23N 


P 2 


T\23N 


-1 


-70.6856984 


-25 


-70.75620346 


-81 


-70.92141286 




-70.6856977(39) 




-70.75620352(11) 




-70.92141291(70) 




-70.686011 




-70.756299 




-70.921481 


1 


-70.6798310 


25 


-70.609519049 


81 


-70.446146678 




-70.6798305(21) 




-70.609519051(97) 




-70.446146655(35) 




-70.680106 




-70.609231 




-70.446044 



Table 5: Comparison for small \p 2 \ < (mi + m 2 + m^) 2 and for mi = 10, m 2 = 20, = 
100. In each box the first entry is the multiple series value of ||, the second entry is our 
result (the error on the last digits is in parenthesis) the third entry is from the numerical 

Our value of p 2 corresponds to — p 2 in 



integration in Table 7 of 11 



[@ and to s m 11 



presented for small s, equal to our p 2 , \p 2 \ < (mi + m 2 + m^) 2 and for mi = 10, m2 = 
20, r«3 = 100. They are repeated here in Table |5|, where in each box the first entry comes 
from the multiple series of with a large number of terms, the second entry is the present 
result, the third entry is from the numerical integration of [l~T|. Again we have excellent 
agreement with the multiple series of 0, while the accuracy of the numerical integration 
of |TTJ is within a few ppm, inside the relative 10 -5 precision declared there. 

In Table || we report the results of Table 2 of || for the combination Xi 23 Ar, for large 
and negative p 2 (i.e. \p 2 \ > (mi + m 2 + m^) 2 , the value of p 2 here corresponds to — p 2 in 
H), so that this time there is also an imaginary part. In each box the first entry comes 
from the multiple series of 0, with a large number of terms, the second entry is the 
present result. Again we have excellent agreement with the multiple series of in most 
of the cases, in few cases there is a deviation of two times the assigned error, that we 
attribute to our procedure of approaching the zero mass. In the seventh box we assign an 
error also to the multiple series, because, although each sum is taken up to 70 terms, the 
results are not yet stable. The assigned error is the difference with the sums taken up to 
60 terms. We attribute the difficulty to the chosen value of p 2 = —150, which is too near 
to the threshold value -(3 + 4 + 5) 2 = -144. 

In Table |7| we report the results of Table 2 of for the combination Ti 23 tv, for large 
and positive p 2 (i.e. \p 2 \ > (mi + m 2 + m 3 ) 2 ; our p 2 corresponds to — p 2 in 0). In each 
box the first entry comes from the multiple series of 0, with a large number of terms, 
the second entry is the present result. Also here we have excellent agreement with the 
multiple series of || in most of the cases, in few cases there is a deviation about two times 
the assigned error, that we attribute to our procedure of approaching the zero mass. 
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777 1 


T77o 


'"'3 


T? 1 
P 


Re Tyizn 


Illl Xi23iV 


2 


3 


2 


-80 


0.587432001 

fl 587431 990f41 "1 


-11.262835755 
-11 9fi98S5744f94 N l 


3 


3 


3 


-100 


-1.28284949 

-1 98984943f 1 61 


-20.899600723 
-9D 89960D741 fl f)1 


2 


3 


4 


-100 


-0.3286481685 

-D 3986481 687f1 61 


-11.84587606309 
-11 845876063D9f6nl 


3 


4 


4 


-150 


-1.26795173 

-1 967951 73f91 1 


-26.491194705 
-96 491 1 94fi94f41 1 


2 


4 


3 


-150 


1.5662482672 

1 566948967Df1 41 


-11.9689438774 
-11 9689438778f 1 91 


3 


3 


4 


-150 


0.9865824304 
986582431 2 fl 9") 


-16.10213970663 
-16 10213970687(77") 


3 


4 


5 


-150 


-4.7638745(12) 
-4.7638748416(32) 


-29.601246304(26) 
-29.60124631204(68) 


2 


3 


4 


-200 


1.6960823345 
1.6960823350(17) 


-6.02417248918 
-6.02417248906(81) 


3 


4 


4 


-200 


1.86355967 
1.86355979(34) 


-20.39852237 
-20.39852240(12) 


4 


4 


4 


-250 


2.64395201 
2.64395222(31) 


-27.60904430 
-27.60904439(18) 



Table 6: Comparison for p 2 large (|p 2 | > (mi + m 2 + m^) 2 ) and negative. In each box 
first entry is the multiple series value of Table 2 of 0], the second entry is our result (the 
error on the last digits is in parenthesis). Our p 2 corresponds to — p 2 in ||. 



irti 


m 2 


m 3 


p 2 


T\23N 


mi 


m 2 


m 3 


p 2 


T\23N 


2 


2 


2 


50 


-3.728125558 
-3.728125610(46) 


2 


4 


3 


100 


-7.1836810855 
-7.1836810854(54) 


3 


3 


3 


100 


-8.79126989 
-8.79127000(11) 


3 


4 


3 


120 


-12.430997190 
-12.430997250(30) 


3 


3 


4 


150 


-7.0830520665 
-7.0830520661(28) 


3 


4 


4 


150 


-10.85647158 
-10.85647158(26) 


3 


4 


3 


150 


-11.361931056 
-11.361931121(31) 


3 


4 


4 


200 


-9.64611359 
-9.64611360(29) 


2 


3 


4 


200 


-3.3018636831 
-3.3018636830(21) 


4 


4 


4 


250 


-13.57188440 
-13.57188463(21) 



Table 7: Comparison for p 2 large (\p 2 \ > {mi + m 2 + m 3 ) 2 ) and positive. In each box 
first entry is the multiple series value of Table 2 of 0, the second entry is our result (the 
error on the last digits is in parenthesis). Our p 2 corresponds to — p 2 in 0. 
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5 Conclusions. 



We propose to solve numerically, by means of the Runge-Kutta method extended to the 
complex plane, the system of the differential equations satisfied by the MI related of the 
diagrams, which due to the large number of occurring parameters cannot be calculated 
analytically. 

We apply the method to the study of the simplest non trivial diagram, the general 
massive 2-loop sunrise self-mass, which is already exhibiting a number of intriguing an- 
alytic properties. We obtain, for all the allowed values of the parameters and of the 
external invariant p 2 , very accurate values of the MI within reasonable CPU time, in 
good agreement with the results already present in the literature. 

The method can be naturally extended to the other self-mass diagrams of the same 
order, like the 2-loop 4-propagator self-mass diagram for which the differential equation 
is already known [|TS] . 



The extension to higher order self-mass diagrams will only increase linearly the number 
of the MI and of the differential equations in the system, while the growing of the number 
of parameters is not a problem at all. Also the extension to diagrams with three or more 
external legs, which means multi variable cases, can be easily envisaged. 

The true difficulty of the method is the need of initial conditions for starting the 
numerical solution of the differential equations; clearly the initial conditions have to be 
provided by an independent method. Of special values are, in this respect, the special 
points (such as 0, oo, thresholds and pseudothresholds) where an analytic calculation is 
easier and sometimes possible, as in the case discussed in this paper. When the special 
points are used, also the first derivative of the MI have to be provided as an independent 
input to the Runge-Kutta approach, but that is analytically a relatively simpler task, 
amounting to an iteration of the expansions provided by the differential equations. 

Acknowledgments. We thank Sandro Rambaldi for his invaluable advise on the use 
of Runge-Kutta method to solve differential equations. 

One of us (HC) is grateful to the Bologna Section of INFN and to the Department of 
Physics of the Bologna University for support and kind hospitality. 
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A Corrections to some analytic formulae of [|5|, @] 



We report here for completeness the correct form of the formulae, which are wrongly 
reported in our previous publications []5], K and are used here to obtain the numerical 
values of the MI at pseudo-thresholds and threshold. 

In section 5 of f| there are three misprints. In Eq.(41) the factor ^ in front of the 
integral should be missing. In the first line of Eq.(44) \og(y) should read log(ys) and in 
Eq.(47) there is a missing factor 4 in front of Z 2 . 

In || there is one misprint in Eq.(45): the second line from the end should be of the 
opposite sign ('— ' — +'). In Eq.(31) of ]8| the expression for Z 3 (mi, m 2 , m 3 ) should be 
symmetric in all the masses, so it becomes 



I 3 {m u m 2 , m 3 ) = %{m Xi m 2 , m 3 ) + %{m x , m l , m 2 ) - X 3 (m 2 , mi, mi) , (16) 

with 



X 3 (mi,m 2 ,m 3 ) = v /mTm 2 " / dm 3 



Jm${m\ + m 2 + m 3 ) 



log 






\mij 



log( 



1112 



m 3 + mi m 3 + m 2 



ijlog M* " *i) " M* " *2)] + log (^P) + *a) - M* + *i) 

+ log(t - h) [2 log(l - *i) - log(ti)] - log(t - t 2 ) [2 log(l - t 2 ) - \og{t 2 )\ 
- log(t + h) [2 log(l + h) - log(-ti)] + log(t + t 2 ) [2 log(l + t 2 ) - log(-t 2 )] 

-2 Li 2 (*=±) + 2 Li 2 + Li 2 (-^ - Li 2 ( 1 ~ h 



where the expressions of the logarithms account now properly for their imaginary part. 
Consequently the Eq.(42) of [RJ should be replaced by 



6-- = vrf-l + llog(2)Uici 2 (* 
32 V 32 32 6V V 8 V2 
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